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1 Introduction 



For most nonlinear wave equations arising in physical applications, their solitary wave so- 
lutions can be obtained only numerically. The recent interest of the research community 
in such applications as Bose-Einstein condensation and light propagation in nonlinear pho- 
tonic lattices has led to a number of publications where numerical methods for obtaining 
solitary waves in more than one spatial dimension were studied. Most of these recent stud- 
ies focus on the so-called imaginary-time evolution method (ITEM), also referred to as the 
normalized gradient flow method [H [2l [3l H] . In this method, one seeks a solitary wave with 
a specified power, or L2-norm, by numerically integrating the underlying nonlinear wave 
equation with the evolution variable t being replaced hy it (hence the name 'imaginary- 
time'). A key step of this technique is the normalization of the solution's L2-norm to a given 
value at each iteration; it is this step that ensures both the convergence of the method (un- 
der known conditions [4J) and the fact that the solitary wave so obtained has a specified 
powei0. Since the power P and the propagation constant /x of a solitary wave are related 
by a beforehand unknown dependence P = P{f^), then the propagation constant in the 
ITEM cannot be specified and is instead computed using the available approximation to 
the stationary solution at each iteration. 

In some applications, it is more convenient to seek a solitary wave with a specified 
propagation constant rather than with a specified power. This is the case, for example, 
in nonlinear photonic lattices, where the value of the propagation constant conveniently 
parametrizes the localized solution within a spectral bandgap. One numerical technique 
that can be used in this case is the Newton's method or any of its modifications (see, 
e.g., [5]). While this method is known to be very fast and also to be able to converge 
to both fundamental and excited-state solitary waves, it also has drawbacks. First, when 
applying the Newton's method in more than one spatial dimension, one has to invert a 
matrix which is not tridiagonal. To do so time-efficiently, one needs to use one of the 
alternating direction implicit methods, which require a certain programming effort. Second, 
the Newton's method often uses a finite-difference discretization of the underlying equation, 
in which case the accuracy of the obtained solution is only polynomial in Ax, where Ax 
is the typical step size of the spatial grid. Finally, it has recently been shown that the 
Newton's method may suffer erratic failures due to small denominators [6j. On the other 
hand, the ITEM mentioned in the previous paragraph is free of these drawbacks. Namely, 
the inversion of the matrix representing the differential operator [H 0] is done using the 
Fast Fourier Transform, which is a built-in function in major computing software (such as 
Matlab and Fortran) for one and two spatial dimensions and can be readily extended to 

^ In a modification of the ITEM, proposed in [4], one normalizes the peak amplitude (the Loo-norm) of 
the solitary wave rather than the power. The simulations reported in [4] indicate that this version of the 
ITEM is faster and converges for a larger class of solutions than the original ITEM. 
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three dimensions. Also, since the operator of spatial differentiation is implemented using the 
spectral method, the accuracy of the ITEM is exponential in Ax (provided that the solution 
is smooth). In addition, the ITEM does not have the small denominator issue (although in 
most cases it converges only to a dynamically stable fundamental solitary wave 0]). Thus, 
it would be desirable to have a numerical method that would possess the above advantages 
of the ITEM while allowing the user to compute the solitary wave with a specified value of 
the propagation constant rather than with the specified power. 

Such a method has long been known for a class of nonlinear wave equations whose 
stationary form is 

- Mu + uP = 0, (1.1) 

where u is the real- valued field of the solitary wave, M is a positive definite and self-adjoint 
differential operator with constant coefficients, and p is a constant. For example, the solitary 
wave of the nonlinear Schrodinger equation in D spatial dimensions, 

iUt + V^U + \U\'^U = 0, C/(|x| ^ oo) ^ 0, 

92 Q2 (1.2) 

V" = 

upon the substitution 



= -1^ + • • • + 



satisfies the equation 

which has the form (jl.ip with 



C/(x,t) = e*'^*u(x) (1.3) 
{n-V'^)u + u^ = 0, (1.4) 
M = n-V^ (1.5) 



and p = 3. Here ^ is the propagation constant of the solitary wave. We now describe 
the idea of the aforementioned method, which was proposed in 1976 by V. Petviashvili [7] 
and has been referred to in the literature by his name. Petviashvili proposed the following 
iteration algorithm: 

where Un is the approximation of the solution at the nth iteration. In scheme (II. 6p . the 
operators and M can be conveniently implemented via the Fourier transform, e.g.: 

m 



where 



(1.7) 
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= , An / /(k)e'''"fik. 
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and ^[M] is the Fourier symbol of operator M. For example, for the operator (jl.Sp . 
T[M] = /i + k^. Also, the inner product in (11. 6p and in what follows is defined in the 
standard way: 

roo 

{f,g)= r(x)5(x)(ix. (1.9) 



(In fact, Petviashvili stated his iteration scheme using the inner product in Fourier space 
rather than its equivalent form (jl.Op .) Note that the quotient in the parentheses of ()1.6p 
equals unity when Un = u, the exact solitary wave; yet the presence of this quotient ensures 
the convergence of the Petviashvili method when the value of the exponent 7 is taken to be 
in a certain range. Namely, in [7j, where he considered the particular case p = 2, Petviashvili 
also formulated a mnemonic rule which yields, for any p, the value 

7 = ^, (1.10) 

p — 1 

for which the fastest convergence of the iterations p.6p occurs. The origin of this optimal 
value of 7 and the convergence conditions of the Petviashvili method for Eq. (jl.ip were 
rigorously established recently in Ref. [8]. 

The Petviashvili method possesses the two advantages of the ITEM which were men- 
tioned in the second paragraph of this Introduction. Namely, the convenience of its im- 
plementation does not depend on the number of spatial dimensions, and its accuracy for 
a smooth solution is exponential in Ax. Moreover, the Petviashvili method, when it con- 
verges, is quite fast. For example, in the case of the one-dimensional nonlinear Schrodinger 
equation (jl.4p . if one starts with the initial condition uq = e~^^ , one reaches the exact 
solution with the accuracy of 10-1° in just over 30 iterations. Here and below we define the 
accuracy as 

En = — " ' ^— ^ . (1.11) 



{Un,Un) 

Recently, a number of studies reported various extensions of the Petviashvili method to 
equations that are of a form different than (II. ip . In Refs. [9j and [lOj, ad hoc modifications 
of the Petviashvili method were proposed for the following stationary wave equations, arising 
in the theory of nonlinear photonic lattices: 

V'^ u + Vo{cos'^ X + cos^ y)u + = fiu , (1-12) 

V U- — :3 ^— :r = /xu . (1-13) 

1 -|- v/o(cos^ X + cos^ y) + 

In Ref. [11], another ad hoc modification of the Petviashvili method was proposed for the 
so-called generalized Gardner equation, which has a mixed quadratic-cubic nonlinearity: 

(l-dl- adl + d~'^dfj u-u^ + bu^ = o > 0. (1.14) 

However, it is not straightforward to generalize the approaches of Refs. [9l [lOl [H] to 
equations with an arbitrary form of nonlinearity. A different, systematic, modification of 
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the Petviashvili method was proposed in Ref. [12]. This method, referred to in [12] as the 
spectral renormalization, can be extended to equations with arbitrary types of nonlinearity 
and also to systems of coupled equations. One can show that for a single equation with power 
nonlinearity (jl.ip . the spectral renormalization method reduces to the following scheme: 



which is slightly different from the original Petviashvili method (jl.6p . (|1.1U|) . (The original 
Petviashvili form (jl.6p of the spectral renormalization method can be restored if one makes 
a simple modification in Eq. (6) of Ref. [12j.) Moreover, it can be verified that for 
equations that contain a power-law nonlinear term and a potential, as, e.g., Eq. (|1.12p . 
the spectral renormalization method with the slight modification mentioned in parentheses 
above reduces to the method of Ref. [9]; see Example 3.2 in Section 3.2 for more details. 
However, it is not known under what conditions the spectral renormalization method, as 
well as the aforementioned methods of Refs. |9l \TU[ lllj . would converge for a general 
equation or a system of equations. Also, as a minor computational issue about the spectral 
renormalization method, we note that it would require some nontrivial programming effort 
to apply it to equations with a non-algebraic nonlinearity, e.g., to 



In this paper, we present a generalization of the Petviashvili method which can be 
applied to a wide class of nonlinear wave equations (including, e.g., (|1.16p l to obtain some 
of their solitary wave solutions. The idea of this generalization is based on the analysis of 
the original Petviashvili method found in Ref. [8j. We also show how our method can be 
applied to systems of coupled nonlinear equations. The only restriction on the underlying 
physical problem is that it be Hamiltonian. The approximate convergence conditions of our 
method for a single equation are stated and discussed, and they can be straightforwardly 
generalized to the case of several coupled equations. 

The main part of this paper is organized as follows. In Section 2, we first recast the 
original Petviashvili method into an equivalent form. All subsequent analysis will be carried 
out for that equivalent formulation of the Petviashvili method. Also in Section 2, we give 
a summary of the results of Ref. [8] in the form that will be suitable for a subsequent 
generalization. This generalization for a single wave equation is presented in Section 3. 
There, we also give examples of the applications of the new method. Next, in Section 4, 
we extend this method to systems of coupled nonlinear wave equations and present the 
corresponding examples. Thus, Sections 3 and 4 contain the two main results of this study, 
which we summarize in the concluding Section 5. The paper also contains four Appendices, 
whose purposes are described in Sections 3 and 4. 

The reader who is only interested in the main ideas of the generalized Petviashvili 
method, but not in technical details of its practical implementation, can skip the following 




(1.15) 



V^u -|- sinh u = fiu . 



(1.16) 
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material without sacrificing the understanding: all Remarks in Section 3.1 except Remark 
3.1; the entire Section 3.2; all Remarks in Section 4 except Remark 4.1; the entire Section 
4.2; and the Appendices. 



2 Review of the analysis of the Petviashvih method for equa- 
tions with power-law nonlinearity 

In this Section, we first reformulate the Petviashvili method into a different, yet equivalent, 
form. The precise meaning of the word "equivalent" will be stated shortly. Then we review 
the results of Ref. [8] concerning the convergence of the Petviashvili method for equations of 
the form (jl.ip . The way in which we present these results is different from the way they were 
originally presented in [8]. This reformulation of both the original Petviashvili method and 
the results of [5] will prepare the ground for our generalization of the Petviashvili method 
in Section 3. 

To recast the original algorithm (|1.6p into a different form, let us begin by introducing 
a notation. Denote the stationary equation whose solitary wave we want to find by 

Lon = 0. (2.1) 

Thus, in the case of Eq. (jl.ip . operator 

Lo = -M + vP-^ . (2.2) 

Here operator M has the properties listed after Eq. (jl.ip . and u is the exact solitary wave. 
Let us rewrite the iteration algorithm (jl.6p in the form: 



Un+l - Un 



(.. + M-M-M.. + <])(l+ ^ K,Mn„.) ) 
= (un + M-\L,u)r) fl+ ^;"'^y"M "'-n„, (2.3) 

V / V {Un,MUn} J 

where 

{Lou)n = -MUn + < . 

Next, let us linearize the above equation near the exact solution u by substituting into it 

Un = U + Un, |t6„| <C |n| (2.4) 

and neglecting all terms of order O(n^) and higher. Using the equation, (|2.ip . satisfied by 
the solitary wave n, one obtains the linearized algorithm (j2.3p : 

Un+i -Un= (^r^Lun - 7 j^' Ar , (2.5) 

Ar = l, (2.6) 
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where L is the operator of the hnearized Eq. ()2.ip : 

Lun = {-M + pvF~^)un. (2.7) 

We now interpret Eq. (j2.5p as the exphcit Euler discretization of the following continuous 
linear flow: 

Ur = M-'Lu-^p^u, (2.8) 
{u, Mu) 

where r is the auxiliary (nonphysical) "time" variable. In the last and key step of this 
derivation, we "de-linearize" the above continuous flow: 

{u, Mu) 

where the notation u simply signifies that this variable is the "current" approximation to 
the exact solitary wave u. That is, if one linearizes Eq. (j2.9p via a continuous analogue of 
()2.4p . one will obtain Eq. (|2.8p . Finally, we discretize Eq. ()2.9p in time using the explicit 
Euler method: 

V {Un,MUn) J 

where (Lou)„ is defined after Eq. (|2.3p . Algorithm (|2.1Up with At = 1 is equivalent 
to the original Petviashvili algorithm (jl.6p in the sense that the linearizations of both 
algorithms yield the same pair of equations (|2.5p and (j2.6p . In the remainder of this paper 
we will, therefore, refer to algorithm (j2.10p also as the Petviashvili method. Moreover, the 
generalized Petviashvili methodproposed in Section 3, which is one of two main results of 
this paper, will be based on this reformulated version of the original algorithm (jl.6p . 

Let us point out two reasons why Eq. ()2.10p is preferred over Eq. (jl.6p for the subsequent 
generalization of the method. First, the ability to select the value of the new parameter 
At makes the convergence conditions of scheme ()2.10p more relaxed than those of the 
original scheme (jl.6p . as shown by Eq. (j2.24p below; see also a related discussion about the 
ITEM in [1]. Second, when trying to generalize algorithm (jl.6p . one may encounter the 
situation where the quotient in the parentheses is negative and hence cannot be raised to 
a non-integer power (without making Un+i complex- valued); see, e.g., [HITO]. In contrast, 
algorithm (|2.10p is free of this difficulty. 

We now come to the second part of this Section where we will review those of the 
calculations of Ref. [8] which are essential for our own analysis given in the next Section. 
Specifically, we will exhibit the conditions under which iterations (j2.10p converge to the 
solitary wave u, that is, when the error Un tends to zero as n ^ oo. To find out when this 
occurs, one substitutes the following decomposition of Un into (12. 5p : 

tt„(x) = ann(x) z„(x) , (2.11) 

where a„ is a scalar (i.e., not a function of x) and Zn{x.) is chosen to be orthogonal to Mu 
at every iteration: 

{zn,Mu)=0, or {Mzn,u) = for all n. (2.12) 
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The second orthogonality relation follows from the first one because, by our assumption, M 
is a self-adjoint operator. 

Before we proceed, let us first point out a relation that will be of crucial importance 
both for the remainder of this section and for Sections 3 and 4. Namely, for Eq. (jl.ip . we 
use Eqs. ([231), <^M, and dSI) to obtain: 

Lu = {p-l)Mu, (2.13) 

or, equivalently, 

M-'^Lu= {p-l)u. (2.14) 

Thus, u is an eigenfunction of operator M^^L, which is closely related to the operator on 
the r.h.s. of Eq. (j2.5p . Equation (j2.14p is the key relation mentioned above; establishing 
its counterpart for a more general Eq. (j3.ip below will correspondingly be one of the key 
steps in the generalization of the Petviashvili method in Section 3. Let us note that from 
Eqs. (|2.12p and (|2.13p there follow the orthogonality relations 

{zn,Lu)=0, or {Lzn,u) = for all n. (2-15) 

Here we have used the fact that L is self-adjoint. 

We now continue with the analysis of the evolution of the error n„ with n. Substitut- 
ing decomposition (|2.1ip into Eq. (j2.5p and using relation (|2.13p and the second of the 
orthogonality conditions (j2.15p . one obtains: 

(on+i - an)u + (z„+i - Zn) = M^^Lzn At + a„n(p - 1)(1 - 7) At . (2.16) 

Taking the inner product of this equation with Mu and using the orthogonality conditions 
([2T2]) and ([2T3]) . one gets 

an+i = a„(l + (p-l)(l-7)AT) . (2.17) 

Thus, when 

a„+i = 0, i.e. the component of the error Un+i "along" the eigenfunction u is zero (in 
the order 0{un)), no matter what this component was at the nth iteration. Note that for 
At = 1, formula (|2.18p yields the optimal value (jl.lOp of 7 found empirically by Petviashvili. 

When a„ and a„+i are related by expression (|2.17p (for any 7), the component Zn of 
the error satisfies: 

Zn+i = (l + ATM-^L)zn. (2.19) 



Since L is self-adjoint and M both positive definite and self-adjoint, eigenfunctions ipj of 
M~^L, satisfying eigen-relations 

M-^Lipj = Xjipj , (2.20) 
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form a complete set in the space of square-integrable functions and are mutually orthogonal 
to each other with weight M: 

{ij,,Mi;k) = Sjk- (2.21) 
Then Zn and ^n+i can be expanded over this set: 

Zn= ^J^ni'j, (2.22) 

where Zj^n are the expansion coefficients. The term with ipj = u (see (j2.14p ) is excluded 
from the above sum because z„ is orthogonal to u with weight M (see (|2.12p and (|2.21|) ). 
Thus, if the eigenvalue {p — 1), corresponding to the eigenfunction u, is the only positive 
eigenvalue of M~^L, then z„ is expandable only over the eigenfunctions with nonpositive 
eigenvalues Xj, and the expansion coefficients satisfy 

Z,,n+i = (1 + Ar)Z,- „ . (2.23) 

As long as At is taken sufficiently small to ensure that 

l + AminAr>-l, (2.24) 

then \Zj^n\ as n — > oo, and hence limn^oo l-^nl — 0. Given that — (in the order 
0{un)) at every iteration when 7 is chosen according to (|2.18p . decomposition (j2.1ip implies 
that \un\ — > as n ^ 00. That is, the Petviashvili method, under the above conditions, 
converges to the solitary wave u. 

3 The generalization of the PetviashviU method for a single 
nonlinear wave equation with a general form of nonlinear- 
ity 

This section contains the first main result of this study. Namely, we will show how the 
Petviashvili method can be generalized for an equation of the form 

Lqu = -Mu + F{^,u) = u(|x|^ 00)^0, (3.1) 

where F{x,u) is any real-valued function. In Section 3.1, we will derive and discuss the 
algorithm of this method and in Section 3.2 will illustrate it with examples. 

3.1 Derivation of the generalized Petviashvili method 

Recall that one of the key results of Section 2 was Eq. (j2.13p . It was that relation on 
which the usefulness of decomposition (j2.1ip was based; see the derivations of Eqs. (j2.16p 
and (|2.17p . Therefore, we will seek to obtain a counterpart of (j2.13p for Eq. (j3.ip . To 
make the main problem of obtaining such a counterpart clearer, consider a particular case 
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of that equation that arises, e.g., in the theories of Bose-Einstein condensation and hght 
propagation in nonhnear photonic lattices: 

Lou= -Mu + V{^)u + u^ = 0, m(|x| ^ cx)) ^ 0. (3.2) 

In the aforementioned physical applications, M is given by Eq. (jl.Sp . and V{x) is some 
potential. The linearized operator L in this case is 

L = -M + y(x) + 3n^ (3.3) 

and hence 

Lu = -Mu + y(x)n + 2>u^ = 2u^ = 2(M - ^(x))^ / const • Mu . (3.4) 

Thus, an exact counterpart of Eq. (|2.13p for a general stationary wave equation p.ip cannot 
be obtained. 

As a solution to the above problem, we propose to seek such a positive definite and 
self-adjoint operator N that the counterpart of (|2.13p would hold approximately: 

Lu aNu . (3.5) 

Here both and the constant a remain to be determined. Given such and a, we then 
construct the following counterpart of algorithm (j2.10p : 

n„+i -Un= [N {Lou)n - 7 — — ^ Un] At, (3.6) 

V {Un,NUn) ) 

where 

7 = 1 + ^. (3.7) 

The algorithm given by the iteration scheme (j3.6p . (j3.7p is the main result of this section. 
We will refer to it as the generalized Petviashvili method. All the steps of the analysis in 
the second part of Section 2 can now be repeated, leading to the following approximate (see 
below) convergence condition for this new method: 

Ij operator N~^L has only one positive eigenvalue (which approximately equals a) and if 
the step size At satisfies inequality 1^2.24^ , where now Amin is the most negative eigenvalue of 
N^^L, then the generalized Petviashvili method converges to the exact solitary wave ii(x). 

Two comments are in order here. First, the component of the error Un "aligned along" 
the eigenfunction of N~^L corresponding to the (only) positive eigenvalue will be annihilated 
in the generalized Petviashvili method not completely, as in the original method (jl.6p , (jl.lOp , 
but approximately. This is due to the fact that u and a are no longer the exact eigenfunction 
and eigenvalue of N^'^L, and hence taking 7 according to Eq. (j3.7p does not make in 
()2.1ip exactly zero at each iteration. That, however, is not really required for convergence: 
It is sufficient that la^+il < |a„| for all n, which is a much more relaxed condition than 
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ttn+i = 0; see also Remark 3.3 below. Second, the reason why the convergence condition 
stated above in italic is approximate rather than exact, as the similar condition for the 
original Petviashvili method stated after Eq. (j2.24p . is the following. Since Zn in (j2.1ip is 
no longer exactly orthogonal to Nu, then the corresponding counterpart of Eq. (j2.19|) for 
method (|3.6p . (|3.7p will hold only approximately. Therefore, in principle it is conceivable 
that if the exact eigenvalue A2 of N^'^L is close to zero and negative, the corresponding 
eigenvalue of the linearized operator on the r.h.s. of (j3.6p will be slightly positive (or vice 
versa). However, such cases are expected to be rare in applications of this method. In fact, 
we did not encounter them in any of the equations to which we applied algorithm ()3.6p . 
(I32D. 

We will now show how the operator N and constant a in (|3.5p can be determined in 
an efficient way. It should be noted that we cannot give the most general recipe in this 
regard, simply because there are infinitely many possibilities here, as it will become clear 
as we proceed. Instead, we will consider in detail only one typical case that arises in many 
applications and will show how N can be found for it. At the end of this subsection we will 
also briefly comment on another example of finding A^. 

Suppose that M in Eq. (13. ip is given by (II. 5[) . The simplest ansatz for is then 

N = c-V^, (3.8) 

where c is to be determined from the condition that "vector" Nun be "aligned along" 
"vector" Lun as closely as possible. Therefore, we require that 

{NUn, NUn){LUn, LUn) 

Differentiating the l.h.s. of the above condition with respect to c and setting the result to 
zero, one obtains 

{NcUn,LUn) _ {NUn,LUn) , . 

{NcUn,NUn) ~ (A'U„,An„)' ^ ■ ^ 

where Nc = dN/dc = 1. The substitution of expression (13. 8p into p.lOp yields the value 
for c at the nth iteration: 

^ _ {Un,LUn){V'^Un,V'^Un) - jV"^ Un , LUn) {u„ , V'^ Un) ,^ 

It is straightforward to verify that for equations with power-law nonlinearity (jl.ip with M 
of the form (jl.Sp . Eq. (j3.1ip yields c = fi and hence N = M. 

Now that A^ has been determined from (13. Sp and (13. lip , the approximate eigenvalue a 
in (j3.5p can be found from 

{Un, LUn) , ^ 

{Un,NUn) 

Thus, Eqs. (j3.1ip . (j3.12p . (j3.6p . and (j3.7p provide all the necessary information for the 
implementation of the generalized Petviashvili method. 
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Before commenting on another example of finding N, we will make several remarks 
regarding implementation of Eqs. (13.111) and (13.121) in a code. As noted at the end of 
the Introduction, the reader who is not interested in such technical details may read only 
Remark 3.1 and then proceed directly to Section 4. 

Remark 3.1 There is no apriori guarantee that the constant c obtained from (j3.1ip 
will be positive, as is required in order to make operator N positive definite. However, in 
all of the examples considered below we monitored c„ and observed it being positive as long 
as we started with a "reasonable" initial condition uq. 

Remark 3.2 This concerns the calculation of quantity Lu„ in Eq. (j3.1ip . Note that 
for any number k, 

{Lu)n + K{Lou)n = (-Mn„ + ^^^(x, M„)ii„) + K (-Mn„ + F(x, n„)) 

= Lu + {Lu)uUn + KLUn + 0{u^) 

= Lu + 0{un), (3.13) 

i.e., in the leading order this expression is independent of k. In (I3.13p . = dF/du, 
{Lu)u = d{Lu)/du, and we have used the fact that Lqu = 0. However, in practice, the 
initial condition uq may not be "sufficiently close" to the exact solution u. This will make 
the 0(nn )-correction comparable in size with the first term on the r.h.s. of (I3.13p . which will 
affect the values of c^, a^, and 7^ in ()3.1ip . (j3.12p . and (|3.7p . This, in turn, may prevent 
the algorithm from converging. In our simulations, we found that this can indeed occur. 
Then we found, empirically, that calculating Lun in ()3.1ip by using ()3.13p with k = —1, 
i.e., 

LUn = {Lu)n - {LQu)n = -F«(x, M„)ti„ - F(x, n„) , (3.14) 

greatly increases the range of the initial conditions uq for which the above algorithm con- 
verges. For example, in the case of Eq. ()3.2p . Lun = 2n^. 

Remark 3.3 The approximate eigenvalue a can be calculated by any formula that 
is equivalent to (j3.12p had (j3.5p held exactly rather than approximately. For example, an 
alternative to (j3.12p may be taken as 

{NUn,LUn) .„^^. 

an = r— ^ . (3.15) 

{NUn,NUn) 

However, in all the examples considered below, we found that the convergence rate was the 
same no matter whether (I3.12p or ()3.15p had been used. This is so because the value of a 
affects only the value of 7, which, in its turn, determines how far the ratio (a„-(_i/a„) is from 
zero. But if this ratio is, say, 0.2 instead of 0.05 (or vice versa) due to a slight variation in 
7, this does not affect the convergence rate, since the latter is determined, in most if not 



all cases, by the much slower decay of ||z„|| = y/ {zn, Zn)- 

Remark 3.4 For some equations, a can be quite small (say, on the order of 0.01 or 
less). We encountered such cases among the examples considered in Section 4. A small a 
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yields a large value of 7; see (j3.7p . We found it to be beneficial to artificailly limit such 
large 7's by using, e.g., 

^ = y2 ' (3.16) 

where 7aux is defined by the r.h.s. of (j3.7p and 7max is some large number specific to the 
problem at hand. The reason why such a limiting may be needed is as follows. Since u is 
not an exact eigenfunction of N~^L, is can be represented as 

u = Uiipi + Ujipj , 
i>2 

where ipj are the true eigenfunctions of N^'^L, and the expansion coefficients Uj are such 
that \Uj\ <^ \Ui\ for j > 2. That is, u contains "small pieces" of eigenfunctions other than 
ipi (the latter is the eigenfunction that u approximates). While the value of 7 as given 
by (j3.7|) is chosen so as to annihilate the ^/^i-component of the error it„, it is not intended 
to annihilate any of the other ^j-components. On the contrary, if 7 is "too large", this 
may amplify some of those components, which would then result in the divergence of the 
iterations. Obviously, this could not have occurred in the case of Eq. (jl.ip . since there 
u = ipi and thus all Uj = for j > 2. 

Remark 3.5 Finally, we note that the computation of c„ at every iteration slows 
down the execution of the code because such a computation requires evaluation of the inner 
products {un,Lun) and {V'^Un, Lun), which are not used in the iteration equation (j3.6p 
itself. However, by the same argument as in Remark 3.3, it is sufficient to compute c„, an, 
and 7„ only until the solution reaches some low accuracy (defined by Eq. (jl.lip ). say, 10~^, 
and then carry on the rest of the iterations using the values of c„ and 7„ computed up to 
that moment. 



As a case where a more involved ansatz for N than p.Sp may appear to be more 
appropriate, consider Eq. (j3.2p in two spatial dimensions where Af is given by (jl.Sp . i.e., 
is isotropic in the spatial dimensions, but the potential V{x) is essentially anisotropic. In 
this case, one may expect that an ansatz more general than (j3.8p . namely, 

N = c-{bdl + dl), (3.17) 

would allow the approximate equation (|3.5p to hold with a better accuracy, which, in turn, 
may result in faster convergence of the iterations. However, this turns out not to be so in 
general. Specifically, we used ansatz (|3.17p for finding solitary waves of equation 

(^sech^(3x) — ij u + u'^ = fiu (3.18) 

(where the potential depends on x but not on y) and observed that not only does using this 
more involved ansatz require more coding effort than when using the simpler ansatz (j3.8p . 
but it also leads to slower convergence of the iterations. To conserve the printed space and 
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the reader's time, we do not show an analysis of this case, since it apparently has little 
practical usefulness. 

We now present three examples that demonstrate the validity of our algorithm (|3.6p - 



(|3in) . (iroi) . 

3.2 Examples of application of the generalized Petviashvili method to a 
single nonlinear wave equation 

Example 3.1 We apply our method to equation 



with A = 1 and W = 1. We take the free parameter At = 1 and compute c and 7 until 
the solution reaches the accuracy of 10~^ (see Remark 3.5). The latest computed values 
are c = 1.20 and 7 = 3.71. Then the iterations are continued until the solution reaches the 
accuracy of 10^^*^. The final solution is shown in Fig. [H and a short Matlab code that can 
be used to obtain it is given in Appendix 1. 

The total number of iterations taken by the generalized Petviashvili method is about 
180. (Here and below we quote the number of iterations rounded to the nearest ten. The 
reason is that this number may slightly depend on the size of the computational domain and 
possibly other technical factors.) For comparison, the optimally accelerated ITEM (with 
the corresponding power of the solitary wave being P = 3.0) reaches the same solution 
in about 300 iterations; see Example 9.1 in [3]. The modification of the ITEM where one 
seeks a solitary wave with a specified peak amplitude [1] rather than the power converges 
to the same solution in 130 iterations. The ad hoc modification of the original Petviashvili 
method, proposed for Eq. (j3.19p in [9], takes about 420 iterations to converge to the same 
accuracy [3]. It should be noted that the value At = 1, which we used in this Example, 
does not lead to the fastest convergence of algorithm (|3.6p - ()3.8p . For instance, we found 
that for At = 1.3, the convergence rate of our method is nearly the fastest, and the method 
converges in about 140 iterations. The dependence of the convergence rate on the step size 
At is discussed in the companion paper jT3] . 

Example 3.2 In this Example we present a case where two modifications of the original 
Petviashvili method proposed in Refs. [9] and [12] diverge, but the generalized Petviashvili 
method, proposed in this Section, converges. We seek an anti-symmetric solution of the 
following equation with a double-well potential: 



V^n + Vo(cos^ X + cos^ y)u + = fiu 



(3.19) 



with Vo = 3 nd ^ = 3.7. The initial condition for the iterations is 



(3.20) 





(3.21) 
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for = 1.43. (The solution with this value of the propagation constant has the power 
P = dx = 10 and was originally found in ^ by the ITEM.) As the initial condition, 

we take 

uo = 2x6"^'' +ee"'^', (3.22) 

with e being either zero or 0.001. As in Example 3.1, we compute c and 7 only as long 
as the error exceeds 10"'^; then we use these latest computed values for the rest of the 
iterations. We first set e = in (j3.22p and empirically find that Ar = 1.6 results in the 
fastest convergence (in about 40 iterations) of the generalized Petviashvili method (j3.6p - 
(13. Sh : the iteratively computed parameters of the algorithm are in this case: c = 5.04 and 
7 = 0.21. The corresponding solution is shown in Fig. [2j Next, when we introduce a small 
symmetric component into the initial condition (j3.22p by setting e = 0.001, the iterations 
still converge to that solution, although at a lower rate (in about 170 iterations). 

We now apply to Eq. (|3.2ip the modifications of the original Petviashvili method pro- 
posed in Refs. [9j and [12]. The former of these methods has the form: 

Un+i = (cZ^^-V{x)un - CZ-^ul) , (3.23) 

where M = fi — d'^ and the factor C„ is chosen so that it equals one when n„ is an exact 
solution of (K2T\\ : 

{Urn U^i) 

The constants 7iin and 7ni in (|3.23p are to be chosen empirically; in the choice 7iin = 0.5 
and 7ni = 1.5 was suggested. The method of Ref. pjj for Eq. p.2ip can be shown to reduce 
to the same form (j3.23p . where now 

_ {Un, (-1 + M-V(x))b„) 

and 7iin = 0.5 and 7ni = 1-5. (Unlike in the ad hoc method of Ref. [9] where these values 
of 7iin and 7ni were "guessed" , in the method of Ref. [12j these values can be derived from 
Eqs. (5) and (6) of that paper.) For the solution of (|3.2ip with ^ = 1.43, both these 
methods converged in about 20 iterations when started at the initial condition ()3.22p with 
e = 0. However, they both diverged for e = 0.001, in contrast to our generalized Petviashvili 
method, which converged for either value of e. 

Example 3.3 In this Example, we show that the calculations of Sections 2 and 3 for 
the optimal value of 7 can be carried out even when the nonlinearity of the equation is 
nonlocal. Consider a stationary wave equation 

l-\yAu = uV-^{^u^\. (3.26) 



2 \ 5x2 
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This equation, which has a nonlocal operator V ^ on the right hand side, can be rewritten 
as a system of two local equations: 

1 



1 — -V^ ] u = uv, 



(3.27) 



This system arises as the small-field approximation in the theory of light propagation in 
photorefractive media; see, e.g., [13]. Let us note that Eqs. ([3.27P cannot be handled by 
the method described in Section 4 below because the corresponding linearized operator is 
not self-adjoint. However, the original nonlocal Eq. ([3.26P can be handled by the method 
of Sections 2 and 3. To that end, let us linearize this equation near the exact solution u: 

Lu^-(l-^vAu + nV-' (-^uA + (^2un\ . (3.28) 



Then, using ([32SD, 



Lu = 2nV-2 ( = 2 - ) n . (3.29) 



This is Eq. ([2TB]) with p = 3 and M = (l - ^V^) ; hence, from ([ITUD . 7opt = 3/2. It is this 
value of 7 which the authors of [H] used (without justification) in the original Petviashvili 
algorithm applied to system ([3.27p . 

4 Generalization of the Petviashvili method for coupled non- 
linear wave equations 

Here we will first show how the generalized Petviashvili method of Section 3 can be ex- 
tended to obtain solitary waves in Hamiltonian systems of coupled nonlinear equations. 
Then we will present the corresponding examples. To make the essential details of our 
technique clearer, we will focus on the case of two coupled equations, while commenting on 
the extension to three and more equations in Appendix 3. 

4.1 Derivation of the generalized Petviashvili method for coupled equa- 
tions 

Consider the following system of equations for the real-valued components u and v of the 
solitary wave: 

V M22 ) \ V ) \ F2(X,U,W) ) \V ) \0 J \^\^oo\v J \0 

(4.1) 

where Mn and M22 are self-adjoint positive definite operators. (Whenever symmetric (see 
below) off-diagonal terms M12V and M12U are present in the first and second equations. 
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they can always be removed by a hnear transformation.) The restrictions on functions Fi^2 
become clear when one considers the linearized operator, L, of Eq. (14. Ih : 



V I \ F2,u -M22 + F2 . I \ V 



(4.2) 



where Fi^u = dFi/du, etc. Recall that the linearized operator L played the key role in the 
analysis of Section 2; in particular, it was crucial for the derivation of Eqs. (j2.15p - (|2.17p 
and the discusion following Eq. ()2.19p that L was self-adjoint. Similarly, to carry out that 
analysis for the coupled Eqs. ()4.ip . we require that L in (j4.2p be self-adjoint. This yields 
the condition 

Fi,^ = F2,u ■ (4.3) 

Thus, our method will be applicable to systems of form (j4.ip where Fi and F2 satisfy 
condition (j4.3p . Note that these functions may contain nonlocal operators as in Example 
3.3 above. 

Our plan now is as follows. We will first present a generalization of the linearized 
continuous flow (j2.8p . then will comment on it, and, finally, will state the vector counterpart 
of the "delinearized" algorithm (|3.6p . The extension of (|2.8p to the vector case is: 



u 



7fc = l + — 7—, ak = -7^—-— , A; = 1,2, 4.5) 

where N is a self-adjoint, positive definite matrix operator, whose form will be discussed 
shortly, and and are the approximate eigenvectors and eigenvalues of N~^L: 

Lgfc « OfcNefc. (4.6) 

The analysis of convergence of Eq. (j4.4p proceeds along the lines of the corresponding 
analysis in Section 2 with one minor modification: In the derivation of Eq. (j4.5p for 7^, one 
needs to use the (approximate) orthogonality of ei and 62- 

(ei,Ne2) = 0. (4.7) 

Condition (j4.7p follows from (j4.6p and the fact that both L and N are self-adjoint. 

Now we discuss the computationally efficient choice of operator N and vectors e^. It is 
this choice that makes the generalization of the method of Section 3 to the case of coupled 
equations nontrivial; hence it constitutes an important technical result of this Section. For 
the simplicity of presentation, we assume that both Mn and M22 in (|4.2p have form (jl.Sp . 
with possibly different ^'s. (The extension to a more general form of these operators is 
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straightforward and one instance of it is given in Example 4.3 below.) Then, the form of N 
that we advocate, and which we used in all of the examples presented in Section 4.2, is 

N=f^' M, Nk = Ck-bkV\ k=l,2. (4.8) 

One may wonder if the more general form that includes (symmetric) off-diagonal terms ci2 — 
6i2V^ would result in a more efficient method. The answer, based on our experimentation 
with both this more general form and the simpler form ()4.8p . is negative. First, the coding 
of the part of the computer program that would calculate all of the coefficients Cfc, bk, and 
ci2i bi2 for the more general form of N is considerably more tedious than the corresponding 
coding for the simpler form ()4.8p . This part of the program would be difficult to debug had 
a mistake in it occurred. Moreover, the simplicity of the original Petviashvili method, which 
is one of its main advantages over the Newton's method, would be compromised by this 
coding issue. Second, in our simulations we also found that, in some cases, unless the initial 
condition {uq,vq) is "very" close to the exact solitary wave {u,v), then the N calculated as 
a full matrix may turn out not to be positive definite, which would result in the divergence 
of the iterations. On the other hand, we verified that the simpler, diagonal form (j4.8p does 
not have either of the above drawbacks. 

Let us now show how ci.2 and 61^2 in (|4.8p can be computed while assuming a general 
form of the eigenvector ei, and then will argue that one can and should take ei = {u,v)'^. 
Let 




where each of Lij is a self-adjoint operator. As in Section 3, we require that 

(Nei, Lei)2 

WT = max 4.10 

(Nei, Nei)(Lei, Lei) 

and then find the equations for ci_2 and 61^2 by setting the derivatives of the l.h.s. with 
respect to these parameters to zero. Thus, similarly to (j3.10p . one obtains a system of four 
equations: 

(N^ei, Lei) ^ (Nei, Lei) . 
(N,.ei, Nei) (Nei, Nei)' ' 
5N 

N^ = r = {ci, 02,^1,62} • 

Since the r.h.s. of all these equations is the same, one can obtain a system of three equations 
for the four unknown parameters r by setting the correponding l.h.s. 's equal to each other. 
It is easy to see, by inspection, that this system is linear and homogeneous, and hence it 
produces a solution for {01,02,61,62} that is unique up to multiplication by an arbitrary 
constant. (Such an arbitrariness is expected because operator N is defined by (|4.6p only up 
to an arbitrary factor.) Next, any one of Eqs. (j4.1ip can be taken as the remaining fourth 
equation for {ci, C2, 61, 62}. We verified that such an equation is satisfied identically for the 
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solution {ci, C2, 62} determined from the aforementioned linear homogeneous system of 
three equations. 

The solution of that system can most easily be found as follows. Equating the l.h.s.'s of 
Eqs. (j4.1ip with r = to those with the corresponding r = hk yields: 

(V^Cfci, Ej=i LkjCji) efci - (efci, J2j=i Lkjeji) V^ekij , V^ekiJ 

" — — — 7~"7 \ \ — 5 ^ — I5 2. 



bk 



( (V^Cfci, ^2%! Lkjeji) Cfci - (cfci, Ei=i Lkjeji) V^efei) , e^i 



Then, equating the l.h.s.'s of (j4.1ip with r = ci and r = C2 yields 

^2 _ (en, (ki - V2)eii) (e2i, Ej=i-^2ieji) 



(4.12) 



(4.13) 



h (621, (k2 - V2)e2i) (en, J2j=i Lijeji) 
The pseudocode for the time-efficient computation (i.e., a computation that avoids repeated 
evaluation of the same quantities) is presented in Appendix 2. As in Remark 3.5, we note 
that C1.2, fci,2 and the corresponding values of 01^2 and 71^2 need only be computed up 
to the moment when the iterations approach the exact solution with some relatively low 
accuracy (say, 10~^). The remaining iterations, up to a higher accuracy, can be carried out 
with those latest computed values of these parameters. 

Remark 4.1 Let us reiterate that the above algorithm of finding the coefficients of 
operator N, which can be straightforwardly generalized to any number of coupled equations 
(see Appendix 3), is one of the main results of this Section. The key part here is that a 
unique set of these coefficients can always (except, maybe, in some pathological cases which 
we never encountered) be found by solving a linear system of equations. 

We now discuss the choice of the eigenvectors ei and e2. First, we note that since these 
eigenvectors enter Eq. (j4.4p on equal footing, it might seem that it would be "more correct" 
to replace the l.h.s. of (j4.10p by 

V (4 14) 

(Nefc, Nefc)(Lefc, Le^) 

However, this is not so because, in particular, the corresponding counterpart of (j4.1ip 
becomes a truly nonlinear system for {ci, C2, 62} and hence cannot be easily solved. 
Therefore, we continue to use the results obtained from ()4.10p . Next, a reasonable, although 
not the most general, choice for ei is 



gi = . (4.15) 

V P21V J 

Then 62 is sought in the form 

e-, = ( "''l ) , (4.16) 
where pi2 is determined from the orthogonality condition (14. 7p : 

Pi2 = -P2i 7 r^, (4.17) 

(m, Niu) 
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where A^i and are found from (j4.8p . (j4.12p . and (j4.13p for each given value of p2i- 

The issue is then to determine coefficient p2i- This can be done by imposing the require- 
ment that quantity (j4.14p . which is a nonhnear function of p2i, be maximized with respect 
to that coefficient. It can be shown, with some effort, that this nonhnear optimization prob- 
lem can be solved time-efficiently, i.e. without repeated evaluation of the inner products in 
(j4.14p . We performed several experiments with the Examples reported in Section 4.2 and 
concluded that simply taking 

P21 = 1 (4.18) 

instead of solving the optimization problem for that coefficient was the optimal choice, 
for the following reasons. In many cases, we empirically found that the "optimal" value 
for p2i was close to one, and hence the considerable complexification of the code needed 
to compute that value did not justify the obtained improvement of the convergence rate 
by just a few percent. Moreover, in some examples we found that the iterations were 
initially selecting a value of p2i that was not close to (j4.18p . and then they would quickly 
diverge. (This probably occurred when the initial condition was not sufficiently close to the 
exact solution.) On the other hand, setting p2i according to (14.18(1 always resulted in the 
convergence of the iterations. Thus we conclude that taking the eigenvectors ei^2 according 
to Eqs. ()4.15p - (l4.18p constitutes the optimal practical choice. The results presented in 
Section 4.2 justify the validity of this choice. 

We now state the algorithm of the generalized Petviashvili method for coupled nonlinear 
wave equations, which is obtained by "delinearizing" Eq. (14. 4p : 




(4.19) 

where e^^n are computed using the components Un,Vn at each iteration, and N and are 
computed iteratively until the solution reaches a prescribed accuracy (see Remark 3.5). 
Iteration scheme (j4.19p along with the details of calculation of N and e,fc (Eqs. (j4.8p . (j4.12p . 
(j4.13p . and (|4.15p ~ (j4.18p ) is the main result of this Section. As we noted in the Introduction, 
the reader who is not interested in implementation issues of this algorithm may skip the 
remainder of this Section. 

Remark 4.2 This Remark extends to the case of coupled equations the observation 
stated in Remark 3.2. Namely, to calculate the coefficients of N and the eigenvalues ai_2 
at the (n -|- l)st iteration, one requires the values of Lei^2) where §1^2 are found from (j4.15p 
and ()4.16p using the available values of u„ and f„. Now, the expressions 

LSk + const • Lo ( ^" ) (4.20) 
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are equal to each other up to the order 0{un,Vn) for any value of the constant, and so the 
issue is which of these expressions to use when computing Lei^2- In our simulations, we 
found that computing Le^ at the nth iteration as 

Lek = {Lek)n - (lo (^""^^^ , A; = 1,2 (4.21) 

(i.e. taking in (14.201) const=— 1) results in a sufficiently broad range of initial conditions 
(no,fo) that converge to the solitary wave {u,v). For comparison, taking in (|4.20p const=0 
required the initial conditions to be much closer to the exact solution for the iterations to 
converge. Thus, to compute qi^2 in (|4.5p . we used the expression given by ()4.2ip . Note that 
while for k = 1, this result is an obvious extension of ()3.14p (upon taking into account (j4.15p 
and (j4.18p ). for k = 2 this result is not obvious and was arrived at upon experimentation 
with various values of the constant in ()4.20p . 

Remark 4.3 When system ()4.ip is decoupled, i.e., Fi^^ = F2^u = 0, the approximate 
eigenvalues ai,2 ofN~^L must be equal. Indeed, in this case, from ()4.5p one has: 

{u, Luu) 1 + {{v, L22V)/{U, Liiu)) 



ai 



OL2 



(4.23) 



(u, iVin) 1 + ((v, N2v)l{u, iVm)) ' 
(u, Liiu) Pi2 + ((^, L22v)/{u, Lnu)) 
{u, Niu) ■ pj^ + {{v, N2v)/{u, Niu)) ■ 

Next, using Eq. ()4.13p with L12 = L21 = 0, one obtains: 

{v, N2V) _ 62 {v, {k2 - V^)v) _ {v, L22V) 
{u, Niu) bi {u, (k2 — V^)u) (n, Luu) 

Substituting (j4.23p into (j4.22p . one obtains ai = 02- This fact is, thus, a consequence of 
the coefficients of the entries of N satisfying (j4.13p . 

Remark 4.4 For completeness of this presentation, we note that it is possible to find 
such a form of Eqs. (j4.ip for which the coefficients of operator N, the coefficient pi2, and 
the eigenvalues ai_2 (and hence 71,2) can be obtained analytically (i.e., similarly to how 
the optimal 7 given by (I2.18P is obtained for Eq. (jl.ip ). For the case of two coupled 
equations, we derive the corresponding class of equations in Appendix 4. A particular 
equation from that class is considered in Example 4.3 below. Note that for this class of 
equations, N = diag(Mii, M22), where Mn and M22 are defined in Eq. (14. ip . This is a 
counterpart of the relation N = M for a single equation with power-law nonlinearity, noted 
after Eq. (l3TT]) . 

4.2 Examples of application of the generalized Petviashvili method to 
coupled equations 

The examples presented below are restricted to systems of two coupled stationary wave 
equations. We focused on those examples where the components u and v of the solitary 
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wave have distinctly different amplitudes and widths; this is done to apply as strict as 
possible a test to our method. Also, in all of these examples except Example 4.1, the 
computational domain was a square with the side of Svr and 2^ mesh points along each side. 
In Example 4.1, the side of the square was 127r with 2*^ mesh points per side. 

Example 4.1 Consider a vector generalization of the equation from Example 3.1: 

V^u + 4(cos^ X + cos^ y)u + uiv? + av'^) = /iiu 

(4.24) 

V'^v + 4(cos^ X + cos^ y)v + v{av?' + 'iv'^) = fi2V. 

Here the asymmetry between u and v is provided by two sources: (i) by the different 
coefficients, '1' and '4', in front of the self-nonlinearity terms and, more importantly, (ii) 
by the different propagation constants m and //2- Specifically, we used 

/.ii = 4.95 and fi2 = 6.5. 

In the absence of coupling (u = 0), this corresponds to the solution u being near the edge 
of the zeroth band gap and v being suficiently far away from that edge. Consequently, v 
is significantly "taller" and more localized than u (see, e.g., |l5j). When the coupling is 
present {a > 0), the structure of the composite solution remains qualitatively the same; 
such a solution for 

a = 0.5 (4.25) 
is plotted in Fig. [3l Starting with the initial condition 

where Ai^2 and Wi^2 are listed in Table 1, the iterations (j4.19p with At = 1 |^ takes about 
710 iterations to converge to accuracy of 10"^''. Here and below, the accuracy for two- 
component solitary waves is defined similarly to (|1.11|) : 

f {Un - Un-l,Un - Un-l) , {Vn - Vn-l,Vn - Vn-l)\^^'^ / a 

i^n = -, ^ \ -, ^ • (4.27) 

In all of the examples of this Section, we monitored the following quantities: coefficient 
pi2 (see (gm-dllH])); factors 

(Nefc, Nefc)(Lefc, Lefc)' ^ ' 

which show how close vectors are to the true eigenvectors of N~^L; the eigenvalues a\^2 
(see (j4.6p ); and the coefficients ci, C2, 62 of N (we set 61 = 1 without loss of generality) . These 
quantities are reported in Table 1. In particular, one sees that ei is a closer approximation to 
its corresponding true eigenvector of N^^L than e2 is to its true eigenvector; this is expected 



^This value of At is likely not to be optimal (see Example 3.1). However, our focus here is not to optimize 
the convergence rate but to demonstrate the vahdity of the method. 
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since the coefficients of operator N are computed using ei. Note, liowever, from the reported 
values of Ii,2, that both ei and §2 approximate their respective true eigenvectors quite well. 

To benchmark the performance of the method, we also obtained the solution of the 
uncoupled system, (j4.24p with cj = 0, in two ways. First, we used the vector generalization 
of the Petviashvili method described in Section 4.1. For Ar = 1, the iterations converged 
to accuracy 10^^^ in about 950 iterations. (Let us note, in passing, that the numerically 
found ai^2 agree with Remark 4.6.) As the second way of obtaining the same solutions, we 
solved each of the uncoupled equations (j4.24p using the generalized Petviashvili method for 
a single equation, as described in Section 3.1. The iterations for components u and v took, 
respectively, about 950 and 80 iterations to converge to the accuracy of 10^^*^. Comparing 
this with the number of iterations needed to obtain the solution of the uncoupled system via 
the first method, we conclude that the convergence rate of the vector form of the generalized 
Petviashvili method is determined by such a rate for the more slowly converging component 
of the solitary wave. 

Example 4.2 We now consider a system of linearly coupled nonlinear Schrodinger 
equations: 

V^M + + av = u 

(4.29) 

V^f + + (7U = V . 

This system extends to two dimensions the equations of the so-called nonlinear directional 
coupler [16]. In one spatial dimension, these equations are known to possess symmetric 
(m = v), anti-symmetric {u = —v), and, for a < 0.6, asymmetric (|m| 7^ \v\) solitary waves 
|16j . The smaller the a, the greater the asymmetry between the two components of the 
latter solution. To our knowledge, asymmetric solutions of the two-dimensional system 
(j4.29p have not been reported previously. 

We considered Eqs. (j4.29p with a = 0.5. By trial and error, we found that the initial 
condition (j4.26p with the parameters reported in Table 1 leads the iterations to converge 
to the solution depicted in Fig. HI (It should be noted that this initial condition must be 
quite close to the exact solution in order for the iterations to converge. For example, if one 
takes A2 = 0.4 or A2 = 0.6 instead of 0.5, as in Table 1, then the iterations converge to 
either the symmetric or anti-symmetric solitary wave.) Next, by running the simulations 
and monitoring, at each iteration, the approximate eigenvalues ai,2, we observed that 02 is 
a large negative number (see Table 1). Then, to satisfy the necessary convergence condition 
(j2.24p . one needs to use a rather small step size At. By trial and error, we found that 
Ar = 0.08 results in nearly the fastest convergence of the method for system ()4.29p . 

Note that since 02 < in this example, the iterations would still converge if 72 were set 
to zero. 

Example 4.3 As the last example, we applied method ()A3.2p to a system of equations 
that describe copropagation of the fundamental and second harmonic fields in an optical 
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medium with quadratic nonlinearity: 

V^n + uv = fiiu 

+ = ^^2v, vj = dl + 5dl . 

Multidimensional solutions of this and related systems were considered in quite a few studies; 
see, e.g., a recent paper [l7j and references therein. It should be noted that Eqs. (|4.30p are 
a special case of the system of two coupled wave equations for which all the coefficients in 
the Petviashvili method can be determined analytically (see Eq. (|A4.14j) in Appendix 4). 
In particular, as follows from the last paragraph of Appendix 4, one should have ci^2 = Mi^2 
and 62 = 1- Thus, this example provides a test of whether our method would obtain these 
coefficients correctly, and it indeed did so. Specifically, we took 

5 = 10, ^ll = 1.5, /X2 = 9; (4.31) 

then one can see that the values of ci^2 and b2, reported in Table 1, are indeed as stated 
above. Moreover, the values of pi2 and 01^2 agree with those given by Eq. (|A4.12l) in 
Appendix 4. Finally, we note that the formulae for the calculation of ci^2, ^2, P12, and 
ai^2 are those given in Section 4.1 with one modification: all occurrences of V^f should be 
replaced with V^f . The corresponding solitary wave is shown in Fig. [H 



5 Summary 

In this work, we obtained the following two main results. 

First, in Section 3, we extended the well-known Petviashvili iteration method to find 
solitary wave solutions of a broad class of Hamiltonian nonlinear wave equations with ar- 
bitrary form of nonlinearity and potential function; see Eq. (j3.ip . Our algorithm is given 
by Eqs. (I3.6p - ()3.8p . (I3.11|) . and ()3.12p . The generalized Petviashvili method can be ap- 
plied even when the equation is nonlocal; see Example 3.3. The computational cost of 
this method only slightly exceeds that of the original Petviashvili method, since the (few) 
parameters required to carry out the iterations need to be computed only until the solution 
reaches some relatively low accuracy; see Remark 3.5. 

Second, in Section 4, we extended this method to systems of coupled Hamiltonian wave 
equations. Our main result here was the finding of a way in which all the required parameters 
of the iteration scheme can be computed by explicit expressions, obtained from solving a 
simple linear system of algebraic equations. The algorithm (for two equations) is given by 
Eqs. gUD, gSD, (BH]), and 

Appendices 1 and 2 contain, respectively, a Matlab code illustrating the algorithm of 
Section 3 and a pseudocode for the algorithm of Section 4. Appendix 3 contains an extension 
of the algorithm of Section 4 to three (and more) equations. Finally, Appendix 4 contains 
a collateral result: the form of a system of two coupled equations for which the parameters 
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of our generalized Petviashvili iteration scheme can be found analytically (as in the original 
Petviashvili method for a single equation with power-law nonlinearity). 
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Appendix 1: Matlab code for Example 3.1 

N=2"7; d=10*pi/N; 'h mesh sizes along x and y 

x=[-5*pi:d:5*pi-d] ; y=x; 

[X,Y] =meshgrid(x,y) ; % 2D x- and y-arrays 

kx=2*pi/(10*pi)*[0:N/2-l -N/2:-l]; ky=kx; 
[KX,KY]=meshgrid(kx,ky) ; K2=KX. '2+KY. "2; 
Dt=0.4; % Delta tau 

mu=3.7; % prop, constant of the soliton 

W=3*((cos(X)) .'■2+(cos(Y)) .-2)-mu; % V(x)-mu 

uO=1.5*exp(-(X."2+Y."2)) ; u=uO; 1 initial condition 

norm_Du=l; % initialize E_n defined in (1.11) 

while norm_Du >= 10" (-10) 

u_old=u; fftu=fft2(u) ; ucube=u.''3; DEL_u=real (-if f t2 (K2 . *f f tu) ) ; 

if norm_Du >= 10" (-3) % when E_n > 10" (-3), compute c and gamma 
dVu=2*ucube; u_u=sum(sum(u. "2) ) ; 

u_Lu=sum(sum(dVu. *u) ) ; DELu_DELu=sum(sum(DEL_u. "2) ) ; 
DELu_Lu=sum(sum(dVu. *DEL_u) ) ; u_DELu=sum(sum(u. *DEL_u) ) ; 
c= (u_Lu*DELu_DEL_u-DELu_Lu*u_DELu) / (u_Lu*u_DELu-DELu_Lu*u_u) 
u_Nu=c*u_u-u_DELu; alpha=u_Lu/u_Nu ; gamma=l+l/ (alpha*Dt) ; 
fftNinv=l./(c+K2) ; % Fourier symbol of N"(-l) 

else °/o once E_n < 10" (-3), use previously computed c and gamma. 

u_Nu=sum(sum(c*u. "2-u. *nabla2_u) ) ; 

end 

LOu=DEL_u + W.*u + ucube; 

u=u+Dt*real (if f t2 (f f t2 (LOu) . *f f tNinv) -u*gamma*sum (sum(u . *LOu) ) /u_Nu) ; 
norm_Du=sqrt (sum(sum( (u-u_old) . "2) )*d"2) ; % new E_n 

end 
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Appendix 2: Pseudocode for time-efficient implementation of 
algorithm (14.191) 



Here we suggest an order in which various quantities, required to perform each iteration in 
algorithm ()4.19p . can be computed. Computing these quantities in this order ahows one to 
avoid repeated time-intensive evaluations, e.g., of inner products such as those required in 
(|4.5p . an so on. 

For notational convenience, we denote, in this Appendix only, 

Ui = Un, U2 = Vn, (A2.1) 

where {un,Vn) are the solution's components at the nth iteration. This notation will facili- 
tate the extension of the steps listed below to the case of more than two coupled equations. 
(A minor modification of this algorithm occurring for more than two equations is described 
in Remark 4.5.) In the notations used below, any index (e.g., j) is assumed to take on the 
values from one to the number of equations (two in the case considered in this paper). The 
summation indices (e.g., k in J2k) o'^^i' same range of values. 

The first column of the list(s) below shows which quantity is computed at the given step. 
The second column shows, which equations of the main text and results of which previous 
steps of this list, are used at the given step. 

The first block of step, listed below, is performed at each iteration, irrespective of the 



magnitude of the error. 

V^ufc {MM (A2.2) 

(Lo)jkUk {MM (A2.3) 

Efc(Lo)ifcnfe {MM (A2.4) 

{uj, Efc(Lo)jfenfe) {MM (A2.5) 



The second block of steps, listed below, contains steps that are required for the calculation 
of the parameters of operator N, the eigenvectors e^, and the parameters 7^. These steps 
need to be performed only while the error is greater than a user-defined threshhold (e.g., 

10-3;,- 

see Remark 3.6 and a note after Eq. Ii4.13\ ). 



LjkUk 


il4.91IA2.2l| 


(A2.6) 




{IA2.6I} 


(A2.7) 


{V^Uj, "Ek LjkUk) 


jlA2.2llA2.6lf 


(A2.8) 


Uk), {uk, y'^Uk), (V^Ufc, V^Ufc) 




(A2.9) 


(V^nj, Y.k{^o)jkUk) 


{IA2.21IA2.4II 


(A2.10) 


(^m(^o)jmUm, J2k{^o)jkUk) 


{|A2.4| 


(A2.11) 
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{Lj^Un,, LjkUk) ^^MM (A2.12) 

{Lj„,u^, Ek(J^o)jkUk) ^^MMIMM (A2.13) 

{uj,j:k Ljkeki) ^rT2nrA23i m7l } (A2.14) 

{V\,EkLjkeki) ^[IMIMMIMM (A2.15) 

Kk ^ m^rM:9llA2. 1411 A2. 151} (A2.16) 

Cfc (take 61 = 1) ^ HTTBirmirM:!! A2.14[rMl6] } (A2.17) 

{uk, NkUk) !^{MM\MM\MM (A2.18) 

P12 ^lim rmsi f (A2.19) 

(Nefc, Ngfc) ^ ]n5iri:T6irTOllA2.18llA2.19l ^ (A2.20) 

(Lefc, Lefc) j lCTll A2.121I A2.13II A2.19I} (A2.21) 

{uj,Y.kLjkek2) i^M\IM{MM{M3\{MM (A2.22) 

{v^uj,Y.kLjkek2) ^imirmmsi mToirmgi ^ (A2.23) 

(Nefc, Lefc) |IA2.14irMl5ll A2.17[rA2l9ll A2.22irA2:23] | (A2.24) 

(gfc, Nefc) ^ IMlSlfACT ^ (A2.25) 

(gfc, Lefc) |IA2.14irA2l9llA2.22l} (A2.26) 

afc, 7fc UMl Dossiblv lXTell A2.251 1 A2.26I} (A2.27) 

The last block of steps is again performed at each iteration, irrespective of the magnitude 

of the error. Note that the latest computed results from the second block are used in this 

one, whereever they are required. 

(ej, Efc(Lo),fenfc) mM[MM (A2.28) 

tifc at next iteration ^ ITOl [A^ I A2.17I A2.191 1 A2.251 rA2:27l I A2.28I} (A2.29) 



Appendix 3: Extension of the algorithm of Section 4.1 to any 
number of coupled equations 

For simplicity, we present the details for the case of three equations; for more equations, 
this treatment can be extended straightforwardly. The counterparts of Eqs. ()4.7p and (j4.8p 
for three equations are, respectively: 

(e,-,Nefc) = 0, j, A; = 1,2, 3, j^k gU) 

and 

N = diag(iVi,iV2,Af3), Nk = Ck-bkV^ k = 1,2,3. gS) 

Then, Eqs. (I4.12[) with k = 1,2,3 are unchanged, and Eq. (j4.13p is replaced with analogous 
expressions for bk/bi where index "2" in ()4.13p is replaced with k = 2,3. Finally, Eqs. 
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(|4.15p and (|4.16p are replaced by 



ei 



u 



( \ 



62 



PI3U 



63 



P23V 
W 



V P32W J 

where w is the third component of the sohtary wave. Since the orthogonahty conditions 
(j4.7l ) yield only three constraints for the four coefficients pjk, we impose an additional 
arbitrary constraint, which we take to be simply 



/O32 = 0. 



(A3.1) 



Then Eqs. (gZI), ^M)-, and ([Mil) yield Eq. (lOTl) (with p2i = 1) for P12 and the following 
system for pi3 and P23: 



Pi3{u,Niu) + P23{V,N2V) = -{w, Nsw) 

Pl3Pl2{u, Niu) + P23{v, N2V) = 0. 

System ()A3.2p can always be solved because pi2 7^ 1. 



(A3.2) 



Appendix 4: Extension of the original Petviashvili method to 
two coupled equations 

Here we will derive the form of two coupled equations for which there exist explicit analytical 
expressions for the coefficients ai^2 etc. (see Remark 4.7). 

Using the analogy with the case of a single equation for which the constant 7 in the 
original Petviashvili method is given by the explicit formula (|2.18|) . we seek the two coupled 
equations in question in the form: 

/ 





(A4.1) 



where Mn and M22 are self-adjoint positive definite operators, as before; pkj and qkj, 
k = 1,2, are some constants; and a^j are linear operators (in particular, they may be 
constants). As we pointed out after Eq. (j4.ip . a more general Hamiltonian system with 
off-diagonal terms M12V and M12U in the matrix above, can be reduced to form ()A4.ip by 
a linear transformation of u and v. The key condition which will allow us to determine the 
relation between the exponents pkj and qkj as well as the parameters of the Petviashvili 
method, is that there is no algebraic relation (such as, e.g., u = const ■ v) between the 
components u and v of the solitary wave. 
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= aiN 




(■) 




V ^ J 







First, we require that the linearized operator L of this equation be self-adjoint (see 
(|4.3p ). which yields that for each j (see the key condition above), either 

qijaij = P2ja2j, (A4.2) 

Pij = P2j-l, (A4.3) 

qij = q2j + l, (A4.4) 

or 

qij = and p2j = 0. (A4.5) 

Next, we require that equation 



(A4.6) 



be satisfied. In view of the key condition, and since ()A4.6P is to be satisfied exactly, it is 
intuitively clear (and can be easily shown) that the only possibility for operator N is: 

■ Mil , , , 

Mi2 = 0, ^ N = I I , (A4.7) 

6M22 

where for the moment constant b is arbitrary. Then ()A4.6P and the key condition yield 

Pi,+qi,-l = ai ^^^^j^^. 
P2j + q2j - 1 = aib 

The counterpart of (|A4.6p for the eigenvector e2 (see ()4.6p and (j4.16p ) yields a similar 
system: 

Pl2{pi,-l)+qi, = a2Pl2 ^^^^^^^ 

Pi2P2j + q2j - 1 = a2b 
Eliminating p2j, q2j, and ai from (jA4.3p . (|A4.4p . and (|A4.8p shows that the following two 
subcases are possible: 

(a) : pij + qij / 1 for ah j, and 6=1; 
(6) : pij + qij = 1 for ah j; 

note that in subcase (b), coefficient b is undetermined. Proceeding with subcase (a), we 
substitute Eqs. (|A43l) and (|A44l) into ([A49]) and obtain: 

Pi2(pii - 1) = a2Pi2 (A4 11) 

puipij + I) + qij - 2 = 02- 
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Since this system is to hold for ah j, with pi2 and 02 being independent of j, one concludes 
that this system can be satisfied only for one set of values {pij,qij}, which we therefore 
redenote as {p, q}. Solving then Eqs. ()A4.1ip for pi2 and 02 and Eqs. ()A4.8P for ai yields: 



P12 



p + l' 



ai 



p + q-1, 



a2 



-2. 



Also, from (|A4.2p - (|A4.4p and ()A4.12p one has: 



a2 = -Pi2 ai 



(A4.12) 



(A4.13) 



Now, using the last equation, one verifies that the orthogonality condition (14. 7p is satisfied 
in this case. Thus, the system of two coupled equations for which the parameters pi2 and 
ai^2 can be determined explicitly (by Eqs. ()A4.12p ) is given by: 



Mil 






M22 



V 



(A4.14) 



where a is any linear operator and p, q are constants. 

Similarly, one can show that subcase (b) of ()A4.10p yields the same equation (jA4.14p . 
where q = 1 — p. Setting the value of the free coefficient b to one yields relations (|A4.12p 
and ()A4.13P in this subcase as well. 

Finally, one can straightforwardly verify that the case given by Eqs. (|A4.5P corresponds 
to two uncoupled equations of the form (jl.ip . Thus, the only nontrivial case in which 
the parameters of the Petviashvili method for a system of two coupled equations can be 
determined explicitly is given by Eq. (|A4.14p . (As we stated after Eq. (|A4.ip . it is assumed 
that there is no algebraic relation between the components of the soltary wave.) The 
parameters of the method are given by Eqs. (|A4.12p . and N is given by (|A4.7p with 6 = 1; 
that is, N coincides with the linear operator in ()A4.14p . similarly to what occurs in the case 
of a single equation with power-law nonlinearity, originally considered by Petviashvili [7j. 
Note that for Eq. ()A4.14p , it is not actually necessary to use the eigenvector 62 in algorithm 
(j4.19p (i.e., one can set 72 = 0), because 02 < and the corresponding component of the 
error would decay on its own (provided that the step size At satisfies the constraint (|2.24p ). 
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Table 1 Values of the parameters, noted around Eqs. (j4.26p and (j4.28p in the text, 
for Examples 4.1-4.3. The asterisk next to the value of Ar means that this time step is 
close to optimal. The numbers of iterations are rounded to the nearest ten. 



Equation 


h,2 


Pl2 


"1,2 


Cl,2 

h2 


^1,2 


WI2 


Ar 


Number of 
iterations 


(I4.24P, 
cr = 0.5 


0.99, 0.69 


-66.2 


0.136, 0.0231 


1.03, 14.9 
7.57 


0.6, 1.5 


2.0, 0.4 


1.0 


710 


(|4.24|), 

cr = 


0.98, 0.78 


-12.4 


0.0943, 0.0943 


1.52, 21.5 
11.0 


0.8, 1.5 


1.0, 0.4 


1.0 


950 


(|4.29|), 
a = 0.5 


1.00, 0.74 


-1.06- 10-2 


2.08, -10.1 


0.750, 0.500 
0.162 


2.0, 0.5 


0.7, 0.3 


0.08* 


580 


(|4.30|), 
(14.311) 


1.00, 1.00 


-0.500 


1.00, -2.00 


1.50, 9.00 
1.00 


1.0, 1.0 


2.0, 2.0 


0.7* 


90 
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Figure 1: Solution of Eq. (f3lJ|l with Vq = 3 and /x = 3.7. 




Figure 2: Anti-symmetric solutions of Eq. ()3.2ip with /x = 1.43. 
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Figure 3: Solution of Eqs. (I4.24j) with /ii = 4.95, fj.2 = 6.5, and a = 0.5. Note the different 
vertical scales of the u- and u-components. 
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-5 5 

X 

Figure 4: Solution of Eqs. (I4.29P with a = 0.5 along the x-axis (the solution is radially 
symmetric) . 




-5 ^ X Y ► 5 



Figure 5: Solution of Eqs. (|4.30p . (|4.3ip along the x-axis (to the left of the dashed hue) 
and the y-axis (to the right of the dashed line) . 
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